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ABSTRACT 

We report on a global, three-dimensional GRMHD simulation of an accretion 
torus embedded in a large scale vertical magnetic field orbiting a Schwarzschild 
black hole. This simulation investigates how a large scale vertical field evolves 
within a turbulent accretion disk and whether global magnetic field configurations 
suitable for launching jets and winds can develop. We find that a "coronal 
mechanism" of magnetic flux motion, which operates largely outside the disk 
body, dominates global flux evolution. In this mechanism, magnetic stresses 
driven by orbital shear create large-scale half-loops of magnetic held that stretch 
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radially inward and then reconnect, leading to discontinuous jumps in the location 
of magnetic flux. In contrast, little or no flux is brought in directly by accretion 
within the disk itself. The coronal mechanism establishes a dipole magnetic field 
in the evacuated funnel around the orbital axis with a field intensity regulated 
by a combination of the magnetic and gas pressures in the inner disk. These 
results prompt a reevaluation of previous descriptions of magnetic flux motion 
associated with accretion. Local pictures are undercut by the intrinsically global 
character of magnetic flux. Formulations in terms of an "effective viscosity" 
competing with an "effective resistivity" are undermined by the nonlinearity of 
of the magnetic dynamics and the fact that the same turbulence driving mass 
motion (traditionally identified as "viscosity") can alter magnetic topology. 

Subject headings: accretion, accretion discs - relativity - (magnetohydrodynam- 
ics) MHD 



Introduction 



Astrophysical jets are seen in such a wide range of astrophysical systems that their cre- 
ation must not require particularly unusual conditions. A tentative consensus has emerged 
that jets are a natural consequence of accretion, rotation and magnetic fields (for a summary 
of this position, see Livio|[2000 ) . Although the presence of accretion and rotation are a given 
in an accreting system, the presence of magnetic fields has long been considered to be less 
certain. However, since magnetohydrodynamic (MHD) turbulence driven by the magnetoro- 
tational instability (MRI; see Balbus fc Hawley||1991 , 1998) accounts for the internal stresses 
that drive accretion, one is at least assured that the presence of accretion itself implies the 
existence of a tangled magnetic field of some reasonable strength. 

But is that enough? Most jet launching mechanisms that have been developed depend on 
the existence of a large-scale, organized poloidal field, following in general terms the scenarios 



outlined analytically in the influential papers of Blandford & Znajek (1977) and Blandford 



& Payne (1982). In the Blandford-Payne model, a large-scale poloidal magnetic field is 



anchored in and rotates with the disk. If the fieldlines are angled outward sufficiently with 
respect to the disk, there can be a net outward force on the matter. As matter is accelerated 
along the rotating fieldlines, its angular momentum increases still further, accelerating and 
driving an outflow. In this model the power for the resulting jet comes from the rotation of the 
accretion disk. In the Blandford-Znajek model, the jet is powered by the rotating space-time 
of a black hole. (In the case of jets from systems with stars rather than black holes, stellar 
rotation can play a similar role.) Radial magnetic field lines lie along the hole's rotation 
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axis and are anchored in the event horizon. Fieldhne rotation is created by frame-dragging 
and this drives an out-going Poynting flux. Most simulations that have demonstrated jet- 
launching have assumed the existence of a large-scale poloidal (mainly vertical) field in the 



initial conditions (see, e.g., the review by Pudritz et al. 2007). Although such simulations 
have provided considerable evidence that magnetic fields can be effective in powering jets, 
they cannot account for the origin of those fields. 

If, as it seems, some sort of large-scale poloidal field is required for jet production, how 
can such a field be established in the near hole region? One possibility is an accretion disk 
dynamo, presumably working with the turbulent fields generated naturally by the MRI. 
A dynamo is an attractive possibility because it has the potential to be ubiquitous. A 
number of models have been put forward suggesting how some form of inverse cascade might 



occur within the MHD turbulence to generate a large scale field (e.g.. Tout & Pringle 1996 



Uzdensky & Goodman 2008), but these scenarios remain speculative and nothing definitive 



has been seen in global disk simulations done to date. In any case, the strict conservation 
of flux within a volume means that a disk dynamo could produce net flux only through 
interactions with boundaries, in this case either the black hole, or by expelling flux to large 
radius. 

The global advection of net flux is an important process, therefore, whether the disk 
functions as a dynamo or net poloidal field is simply carried in to the near hole region from 
large radius. Even without any dynamo action, a weak net field might become locally strong 
it were concentrated by the accretion process. In the absence of a single field polarity at 
large distance, a near-hole net field could nevertheless be built up due to a random walk 



process as field is accreted (e.g., Thorne et al. 1986). Whether or not net field can be so 
accreted, however, is a matter of some uncertainty. The main concern is that the field would 
not be accreted if, as seems intuitively likely, the field intensity declines outward and the 
rate of diffusion of the field through the matter exceeds the rate at which matter accretes 



(van Ballegooijen 1989; Lubow et al. 1994; Heyvaerts et al. 1996 Ogilvie & Livio 2001). 



This competition between inward advection and outward diffusion is typically described in 
terms of an effective turbulent accretion viscosity Ut determining the accretion rate, and an 
effective magnetic diffusivity, rjt that sets the field diffusion rate. The fate of the field is 
then determined by the effective magnetic Prandtl number Vm = i^t/vt- In this picture, it is 
argued that the relevant scale for effective resistivity is ~ whereas the relevant scale for 
effective viscosity is ~ R; consequently, field can be brought inward only if Vm > R/H. 

It is far from clear, however, that field transport within MHD turbulence can be de- 
scribed in such terms. The picture just described assumes that the fluid motion is inde- 
pendent of the field, yet in Nature the dynamical coupling between fluid and field assures 
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that even weak fields are amplified, and that the effective viscosity Pt is, in fact, the result 
of magnetic stresses. Thus, there is no kinematic regime in which the Lorentz forces can 
be neglected. The resulting MRI-driven turbulence is inherently three-dimensional; indeed, 
three-dimensional studies are essential to generating long-term sustained MHD turbulence 
due to the fundamental restriction of the anti-dynamo theorem. By contrast, the frame- 
work of effective turbulent viscosity and diffusivity is cast in axisymmetry, in the hope that 
suitable time- and space-averaging of the field transport process in a fully turbulent, three- 
dimensional disk can lead to some value of r]t in much the same way that suitable averaging 
of correlations in the magnetic field, normalized to some suitable pressure term (either gas, 
magnetic or total pressure) leads to a value of Vf But the primary mechanism by which 
magnetic fields change the fluid elements to which they are attached is through magnetic 
reconnection driven by the fluid turbulence, not gradual slippage via resistivity. The effec- 
tive viscosity and magnetic diffusivity are derived from making the ansatz that stress and 
reconnection can be modeled as a viscosity or resistivity, i.e., a constant coefficient multiply- 
ing the gradient in orbital rotation rate or magnetic field strength. Whether or not such a 
model is appropriate to accretion flows remains to be determined. Finally, net flux and field 
topology are global concepts; a purely local description of field motion in terms of transport 
coefficients may not be sufficient. 

Numerical simulations offer a promising approach to investigate these and other issues 
related to global field evolution within accretion flows. Of course, simulations entail their 
own set of difficulties: global simulations require adequately fine resolution over a wide range 
of radii in the disk body and throughout a good deal of the surrounding coronal region as 
well. The numerical grid boundaries should be sufficiently remote so as to minimize artificial 
effects on the disk. Although axisymmetric simulations can be of interest, the problem 
ultimately must be done in three dimensions because any axisymmetric calculation fails to 
describe properly essential properties of both the fluid turbulence and the magnetic field 
evolution. Finally, because we need to follow the disk long enough to observe net accretion 
in an approximate inflow equilibrium, the duration of the simulation must be long compared 
to typical local dynamical times in the disk. 

Work along these lines is ongoing. Global, three dimensional black hole accretion sim- 
ulations without initial fields threading the disk were carried out and described in a series 



of papers beginning with De Villiers et al. (2003). In these simulations initial weak dipole 



loops are contained within an isolated gas torus. The subsequent MRI-driven accretion flow 
carries the inner half of the initial field loop (with a single sign for the net vertical field 
component) down to the black hole, creating radial field as the loop is stretched. Differential 
rotation along that radial field leads to rapid growth of the toroidal field, especially within 
the plunging region of the flow. This growing toroidal field causes the ejection of a magnetic 
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tower into the low density funnel region and the establishment of a global dipole field an- 
chored in the black hole. This field can then produce a Poynting flux jet if the black hole is 
rotating. 



Beckwith et al. (2008a) examined several alternative initial field configurations in com- 



parison to the dipolar loop. Quadrupolar initial field loops produce a much weaker funnel 



field. A similar result was found by McKinney & Gammie (2004) for axisymmetric simula- 



tions, and McKinney & Blandford (2009) in a three-dimensional model. A variation proposed 



by Hawley & Krolik (2006) and further elaborated upon by Beckwith et al. (2008a) is a series 



of dipolar loops contained within the accretion flow. Such loops can lead to the creation of a 
temporary net sense of vertical field that can support a significant field in a relativistic outflow 



along the rotation axis for a sizable period of time. Beckwith et al. (2008a) also computed 



a model that began with a purely toroidal field. In this case the MRI creates turbulence, 
field and angular momentum transport, but without creating any overall organization to the 
poloidal field. As a consequence, no funnel field is generated. 

In this paper we continue our study of the influence of field topology on accretion and 
potential jet formation through a three-dimensional general relativistic MHD simulation with 
an initial field configuration consisting of a net vertical field passing through an orbiting ring 
of gas at modestly large radius. We will investigate the evolution of the disk in the presence 
of a net field while determining if such an initial condition can lead to the establishment 
of field configurations conducive to jet formation. From this study we hope to gain more 
insight into whether or not large scale fields can be advected along with the accreting gas 
and to lay the groundwork for further work. 

The rest of this paper is structured as follows: In ^ we describe the methodology and 
the initial conditions for our simulations. In ^ we discuss their results as they pertain to 
magnetic flux motion. Finally, in ^ we summarize our results, compare them with current 
models from the literature, and explain their implications for observations of relativistic jets 
in Nature. 



Numerical Details 



For this work we employ the general relativistic MHD code GRMHD developed in De 



Villiers &; Hawley 


( 


2003 


); 


De Villiers 


( 


2006 



GRMHD has already been used in several 



studies to simulate black hole accretion in three spatial dimensions (De Villiers et al. 2003 



Hirose et al.|[2004l |De Villiers eTar][2005| |Krohk et al.][2005| [Hawley fc Krolikl[2006| [Beckwith 



et al.[[2008a[|b| . GRMHD solves the equations of ideal non-radiative MHD in the static Kerr 
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metric of a rotating black hole using Boyer-Lindquist coordinates. Values are expressed 
in gravitational units {G = M = c = 1) with line element cis^ = Qudt"^ + '^gt<j>dtd(j) + 
grrdr"^ + QeedO'^ + g(f,(f,d(f)'^ and signature (— , +, +, +). The relativistic fluid at each grid zone 
is described by its density p, specific internal energy e, 4-velocity U^^, and isotropic pressure 
P. The relativistic enthalpy ish = l + e + P/p. The pressure is related to p and e via the 
equation of state for an ideal gas, P = pe(T — 1). The magnetic field is described by two sets 
of variables. The first is the constrained transport (CT; Evans & Hawley 1988) magnetic 
field = [ijk]Fjk, where [ijk] is the completely anti-symmetric symbol, and Fjk are the 
spatial components of the electromagnetic field strength tensor. From the CT magnetic 
field components, we derive the magnetic field four-vector, (47r)^/^6^ = *F^'^Ui,, and the 
magnetic field scalar, = b^b^. The electromagnetic component of the stress-energy 

tensor is Tj^^^^ = + U^'U''\\b\\'^ - b'^b". 



2.1. Initial and Boundary Conditions 

In this study we carry out both axisymmetric and fully three dimensional simulations 



around a Schwarzschild black hole. The initial condition, as in previous works (De Villiers 



et al. 


2003 


Hirose et al. 


2004 


2006 


Beckwith et al.||2008a|b). 



field is imposed. We use the torus model described by De Villiers et al. (2003) with an 
adiabatic index of F = 5/3. The initial angular momentum distribution is slightly sub- 
Keplerian, with a specific angular momentum at the inner edge of the torus (located at 



30M) 



U(f)/Ut = 6.10. For this choice of ii^, the pressure maximum is at r ^ 40M 



and the outer edge of the torus is located at r ^ 60M. The angular period at the pressure 
maximum is ~ 250M. With this inner edge radius, the resulting accretion flow should 
have sufficient radial range to permit a detailed study of the advection of vertical flux. The 
choice of angular momentum parameter at the inner edge of the torus, ii^, was made so 
that the initial torus (and hence the steady state accretion flow) would have a scale height 
similar to that of previous simulations (see e.g. Beckwith et al.|[2008b ). The torus is initially 
surrounded by zones filled with zero-velocity gas set to a constant value of density and 
pressure. The ratios of these "vacuum" values to the torus maximum are 2 x 10~^ for the 
pressure, and 5 x 10~^ for the density. Because the code is relativistic, the timestep is always 
limited by the time for light to cross the smallest cell; high Alfven speeds in the low-density 
coronal region do not create any unique difficulty. 

The initial magnetic field is homogeneous and vertical, filling the annular cylinder 
35M < R < 55M (where R is the cylindrical radius). This configuration differs from 
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the vertical field used in some previous studies (see e.g. Koide et al. 1999; McKinney & 



Gammie 2004 De Villiers 2006) which fills all space, not just through the initial torus. By 
isolating the field in this manner, we know that any vertical fiux that ends up crossing the 
black hole event horizon was brought in from field initially passing through the torus at 
large radius, rather than already being present at the horizon in the initial conditions. Also, 
because equilibrium field strengths scale strongly with radius, a global vertical field that is 
strong enough to be interesting at one radius is likely to be overwhelming at some radii while 
negligible at others. We prefer to see the global field strength emerge as an outcome of the 
evolution. In our configuration the initial vertical field is subthermal within the torus and is 
unstable to the MRI. It should be noted, though, that this initial field configuration is not 
in dynamic equilibrium outside of the torus. Also, this is but one very specific realization 
of a torus embedded in large scale field. Future experiments will need to explore a wider 
variety of initial configurations. 

This initial field is computed from the curl of the four- vector potential, i.e., 

FaP = daA/s — d/sAa = Ap^^ — Aafi. (1) 

with only A^ ^ 0. The vector potential corresponding to a uniform (Newtonian) vertical 
field geometry is proportional to the cylindrical radius, i.e., 

A^ (xr sin 0. (2) 

The vector potential for the initial field is 

{Rin R ^ Rin 

rsin^ Rin<R<Rout (3) 
Rout R > Rout 

where Rm and Rout are 35M and 55M. Outside of those cylindrical radii the field is zero. 
Af) was specified so that the ratio of the average thermal to average magnetic pressure inside 
the torus is /3 = 100. At the center of the torus P = 275 and /5 declines smoothly toward 
1 near the edge of the torus. In the surrounding corona, (3 = 0.038. The initial mass, field 
and P distribution is shown in Figure [TJ 

The simulation we will focus on is three-dimensional with 256 x 256 x 64 zones in (r, 6, 0). 
The simulation is labeled "VDOmSd." Given the challenge of adequately resolving such 
models and of doing anything like a convergence study, we also carried out two-dimensional 
simulations at two different resolutions. In two dimensions we used 256 x 256 and 512 x 512 
zones in (r, 6*). The two-dimensional simulations are labeled "VD0m2d" and "VD0h2d" 
where the "m" designates medium resolution and the "h" designates the higher resolution 
case. 
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For all models the radial grid extends from an inner boundary at = 2.104 to an 
outer boundary located at Vout = 500M. The radial grid was graded using a logarithmic 
distribution that concentrates zones near the inner boundary. An outflow condition is applied 
at both the inner and outer radial boundaries. In contrast to the zero net flux simulations 
that we have performed previously, magnetic fields pierce the outer boundary in the initial 
state. This net flux is kept fixed throughout the simulation, but otherwise the field at the 
boundary can evolve. When the fluid velocity normal to the boundary is directed outward we 
perform a simple copy of fluid variables into the ghost zones that lie outside of the boundary. 
In the case where the fluid velocity is inward, the ghost zones are filled with a cold, low 
density, zero momentum fluid at the vacuum value so that only cold, very low density gas 
can enter the computational grid. To minimize the influence of the outer boundary it is 
located as far away from the region of interest as is feasible. Any torques that arise from 
the interaction of the outer boundary with the large scale vertical field should have minimal 
influence on the turbulent disk itself. 

The ^-grid spans the range O.Olvr < 6 < 0.9971. This creates a conical cutout that 
surrounds the coordinate axis. Such a cutout will prevent flow through the coordinate axis, 
but it is computationally advantageous to avoid the coordinate singularity. The size of this 
^-cutout is determined by two considerations. First, light travel times across narrow (/)-zones 
near the axis sets the timestep, so the cutout should not be too small. In this simulation. 
At = 2.53 X 10^'^M. The 6 cutout must not be too large; initial vertical field should not cross 
the 6'-cutout before reaching the outer radial boundary at r = 500M. A reflecting boundary 
condition is used along the ^-cutout. The 6 zones were distributed logarithmically so as to 
concentrate zones around the equator. 

The 0-grid spans the quarter plane, < < 7r/2. The grid is uniform and periodic 
boundary conditions are applied. The use of this restricted angular domain significantly 



reduces the computational requirements of the simulation. Schnittman et al. (2006) exam- 
ined the variance in surface density as a function of (f) and found the characteristic size of 
perturbations to be around 25°. Thus, while some global features may be lost by using a 
restricted domain, the character of the local MRI turbulence is captured. 



2.2. Magnetic Flux 

In analyzing the evolution of the large scale poloidal flux we will make use of the poloidal 



flux function, which corresponds to the ^-component of the magnetic vector potential (van 



Ballegooijen 


1989 


Heyvaerts et al. 


1996 



from the simulation data. In GRMHD, we directly evolve the components of the Faraday 
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tensor, F^j, using the alternative form of the induction equation (for further details see De 



Villiers & Hawley 2003): 



dsF^fs + d^Fps + dpFsa = (4) 

The space-space components of F^j^ are identified with the constrained-transport (CT) mag- 
netic fields via = [ijk]Fjk: 

B'- = F^e; l3' = Frf, B* = Fgr. (5) 

This identification allows us to write the induction equation in the familiar form 

dtB' - dj [V'B^ - B'V^) = 0, (6) 

where = W/U^. The tensor F^^i, is related to the magnetic induction in the fluid frame, 
by the relation 

Ffiu = ^aPpLV -B" f^^, (7) 

where e^^s^ = \[ivb-i\. 

The Faraday tensor can also be written in terms of a vector potential: 

Fctfi = daAj3 — dpAa = Ap^a — Aa^fi- (8) 

It is convenient to define an azimuthally-averaged poloidal flux function, A^. This can 
be derived from the azimuthally-averaged CT field data by noting that the total (spatial) 
derivative of is 

dA^ = ^dr + ^de. (9) 
or o9 

/^From the definition of the Faraday tensor, we have that 

F(j)e = Ae^fj) — A^^g; Fr^ = A^^r — Ar^^. (10) 
After azimuthal averaging, this reduces to 

F(j)0 = —A^^g; Fr(l) = Afp^r- (11) 

The space-space components of F^jy are identified with the CT magnetic fields via B^ = 
[ijk]Fjk so that 

B^ = F^0 = —A^^g; B^ = Frtp = A^^r- (12) 

We therefore have that 

dA^ = Fr^dr - F^gde = B'^dr - B' dO. (13) 
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To find the flux from this differential, we could in principle integrate over any poloidal 
surface. For the radial motion of vertical flux through the disk, the surface of greatest 
significance is the equatorial plane. Because this surface is interrupted by the event horizon, 
we must stretch it over half of the event horizon. We will therefore define a flux function 
that is the sum of two parts: the radial flux through the inner horizon boundary. 



and the vertical flux through the equator. 



|6»=7r/2 



rr 

/ dr'B\r\e = Ti/2). 



(14) 



(15) 



Using these functions, we can determine how the accumulated flux through the top half of 
the event horizon depends on polar angle, '^{0), and how much total vertical flux has been 
brought to within radius r. The full flux function along the total path is 



A{r,d = 7r/2) = ^(7r/2) - $(r). 



(16) 



This corresponds to the net flux piercing the surface covering the top hemisphere of the black 
hole horizon and the equatorial plane out to radius r. The full vector potential component 
A^{r, 6) is obtained by integrating the CT variables throughout the computational domain 
starting from r = and 6' = 0. A{r^9 = 7r/2) is therefore identical to A^{r,6 = 7r/2). 
For the three-dimensional simulation, VDOmSd, and are integrated over over the 0- 
dimension prior to computing A(j^{r,6). Given the periodic nature of the direction this is 
sufficient to compute the total flux piercing the volume bounded by the equatorial plane and 
the black hole horizon. We note that this procedure is consistent with the approach adopted 



in the calculations of van Ballegooijen (1989), Heyvaerts et al. (1996), and Reynolds et al. 



(2006) 



3. Global Flux Transport 



3.1. Overall evolution 



The simulation begins with an isolated torus threaded with a column of vertical field. 
The first part of the simulation is marked by relatively rapid evolution of the coronal field 
and the surface layers of the torus. Figure [2] shows the evolution of the poloidal magnetic flux 
and the density distribution over the time period lOOOM <t< 2500M. This process, along 
with the subsequent evolution of the flow is also shown in an animation included in the online 
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edition of this work. We encourage the reader to refer to this material as it clarifies many of 
the issues that we shall discuss at length in the remainder of the text. As noted above, the 
initial state is not in equilibrium within the corona. There is differential rotation between 
the coronal gas and the torus, and the magnetic pressure within the vertical field column 
is not balanced by pressure within the initial atmosphere. The field expands radially, and 
stresses at the boundary of the disk begin to drive low density gas from the torus outward 
along field lines. In the inner half of the torus, low density surface layers move inward in two 
thin accretion streams above and below the equatorial plane, dragging poloidal field lines 
down to the black hole. This behavior was dubbed an "avalanche flow" when it was observed 



in the axisymmetric MHD accretion torus simulations of Matsumoto et al. (1996), and it is a 



common feature to accretion models that begin with vertical fields threading the disk (e.g.. 



Koide et al. 1999). These surface features develop well before field deeper within the torus 
is significantly amplified by the MRI. The surface layers in the outer half of the torus also 
evolve, but here mass flows radially outward along the field lines, carrying the field along as 
it does so. 

By t = 2000M net flux has begun to build up on the black hole horizon, well before any 
significant mass accretion has begun. This flux arrives there by a process we call the "coronal 



mechanism," which we discuss briefly here and will discuss again in greater detail in § |3.2 
A typical field line participating in the coronal mechanism may be traced from the outer 
boundary at high latitude more or less vertically down to the upper surface of the torus, 
where it bends sharply toward the hole near the boundary of the denser region. It goes in 
some distance before doubling back (forming a "hairpin" shape) and returning to the torus. 
There it passes through the equator, and turns inward to mirror its behavior on the other side 
of the equator. A short time later, the apex of each hairpin moves inward and closer to the 
equator, where, in the early stages of this simulation when there is little mass at small radii, 
it can encounter an oppositely-oriented hairpin from the other hemisphere. Reconnection 
ensues, causing the magnetic flux distribution to change essentially instantaneously. 

The four panels of Fig. [2] show these structures in differing ways. In the first panel, 
an isolated hairpin field line can be clearly seen approaching very close to the horizon. By 
the time of the second panel, the bends in this field line have reconnected, attaching it 
to the event horizon and liberating a pair of closed loops. Also in the second panel, new 
hairpins can be seen forming along the top and bottom surfaces of the disk in the region 
20M < r < 30M. These hairpins evolve further and become embroiled in the beginnings of 
disk turbulence (the channel modes) in the third and fourth panels. 

Magnetic field enters the extremely low-density funnel around the rotation axis by an 
offshoot of the coronal mechanism just described. There is initially no field within the 
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funnel. Once horizontal field lines have been brought close to the black hole, their radial, 
near-horizontal components become unstable to the formation of ballooning half-loops that 
rise upward into the funnel. These loops are the initial source of field for the funnel. Only 
the inner half of the rising loop enters the funnel and, as a result, the field direction changes 
sign abruptly at the funnel's outer edge along the centrifugal barrier. 

The coronal process also occurs at radii exterior to the initial pressure maximum of the 
torus. There it acts to move net flux progressively outward. Above the disk surface, the 
field is stretched out radially, but subsequently bends back down toward the equator, where 
reconnection can occur. This transfers the net flux from field lines going through the torus 
(those lines now close in large loops) to vertical field lines at large radius. 

Thus, in the early phase of the evolution the global topology of the field is rapidly 
rearranged through coronal motions without any significant evolution within the main disk. 
Although these features result from the initial conditions, the process illustrates how (in 
principle) coherent large scale poloidal flux can be rapidly moved over large radial distances. 

Although low density surface layers participate in the initial avalanche flow, the dense 
part of the torus is largely unaffected by these coronal motions. By t = 2500M, however, the 
MRI reaches a nonlinear amplitude throughout the torus and the characteristic structures 
associated with the "channel mode" have become visible within the torus. The next stage of 
evolution is illustrated in Figure [3] which shows snapshots every 500M in time from 3000 to 
4500M. Gas in the torus moves radially both inward and outward in the characteristic chan- 
nels as angular momentum is transported along field lines, which are themselves increasingly 
stretched into long radial filaments. The channel modes are unstable to nonaxisymmetric 
perturbations and rapidly break up, resulting in a turbulent disk without prominent radial 
features. In the axisymmetric simulations, these relatively coherent large-scale radial flows 
continue to dominate throughout the evolution. This artificial persistence of coherent motion 
is a significant limitation of two-dimensional simulations with vertical field. 

Accretion into the black hole does not begin until 1650M, and the mass accretion rate at 
the inner edge of the disk reaches its long-term mean for the first time shortly after 4500M. 
The mass accretion rate then continues to vary by factors of ~ 2 around this mean through 
the remainder of the simulation. The disk achieves an approximate inflow equilibrium within 
r ~ 25M from time lO^M onward. Prior to lO^M, the mass in the inner disk (r < 25M) rises 
steadily, but after that time, it changes by about 10%, and without any prevailing trend. 
The mass accretion rate averaged from 10*^M until the end of the simulation at 2 x lO'^M is 
flat to within 10% for r < 25M. The total mass on the grid at the end of the simulation is 
78% of its initial value; of this, 18.8% has gone into the black hole. 
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The divergence-free nature of the magnetic field means that the total flux is rigorously 
conserved except for losses at the boundaries of the computational domain. To measure 
the evolution of the total flux we integrate along the black hole horizon from the axis to the 
midplane and then out along the equator to the outer boundary. Figure |4] shows this integral 
as a function of time. The figure shows that there is no flux lost through the outer radial 
boundary during the simulation, but considerable net flux is deposited on the horizon at the 
expense of flux through the equator. The majority of the flux brought to the horizon arrives 
by t ~ 5000M, which implies that the net horizon field is created mostly by coronal motions, 
not through disk accretion. During the later stages of more or less steady-state accretion, 
the horizon flux grows much more slowly. Relative to its magnitude at t = 5000M, the 
horizon flux increases by only 28% over the next 5000M and adds only another 11% during 
the last lO^M. By contrast, the total mass accreted onto the black hole grows by a factor 
of 9 from 5000M to 2 x lO^M, with a rate of increase that is approximately constant from 
~ lO^M onward. 

The appearance of the overall flow at the end of the simulation is given in Figure |5] 
which shows the azimuthally-averaged density distribution and the (3 parameter, and in 
Figure |6| which shows the vector potential component A,^ at time t = 2 x lO'^M. There are 
several distinct regions. First, there is a low density, strongly magnetized funnel in which the 
field lines are primarily radial. Second, there is a corona characterized by large-scale field 
lines that loop about creating islands and hairpins, with numerous current sheets separating 
closely-spaced regions of oppositely-directed field. Here there are regions of both relatively 
strong and relatively weak magnetization; /3 ranges from ~ 0.1 to ~ 10. A thin region of 
high (3 fluid running along the funnel wall marks a current sheet between field lines directed 
in opposite radial directions as well as the boundary of the region with strong toroidal field. 
Finally, the main disk body remains weakly magnetized, with (3 ranging from ~ 10-100. 



The disk thickness H/R ranges from 0.1 to 0.15, using the vertical density moment (Noble 



et al. 2008) for the definition of H. As can be seen from Figures [3j and [6| substantial parts 
of the disk body (most of the disk inside r ~ 40M) are actually disconnected from the 
large-scale field, their magnetic connections having been broken by the reconnection events 
that transferred field both to the horizon and to the outside of the torus. In the outer half 
of this region (20M < r < 40M), however, some field lines do pass through the equator and 
out through the disk and into the corona. 

This initial magnetic field and torus configuration was also simulated in axisymmetry. 
While the initial evolutions of the two- and three-dimensional models are similar, once the 
MRI sets in within the torus the subsequent evolutions differ dramatically. In axisymmetry, 
the channel modes dominate throughout the simulation, creating coherent radial field lines 
that extend over large distances. These field structures strongly influence the accretion 
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of both mass and magnetic flux. Figure [7] compares the flux through the horizon, the 
mass accretion rate and the specific angular momentum accreted in simulations VDOmSd 
and VD0m2d. The axisymmetric simulation is characterized by very large fluctuations in 
mass accretion rate. Similar large fluctuations are seen in the value of the magnetic flux 
through the horizon. Note the contrast between the net specific angular momentum jm 
carried into the hole by the accretion flow in the two- and three-dimensional simulations. 
For both plots the dashed line indicates the specific angular momentum of the innermost 
stable circular orbit at 6M. The deep minima in the axisymmetric run correspond to strong 
torques created by radial field extending through the plunging region. These extended radial 
fields are associated with the coherent channel flows seen in axisymmetry, and they have an 
obvious strong effect on mass, angular momentum and magnetic fluxes carried into the black 
hole. The three-dimensional simulation is qualitatively different. In three dimensions, the 
channel modes quickly break up, yielding to genuine turbulence. The accretion rate into 
the hole fluctuates, but not nearly as strongly as in axisymmetry. Similarly, only in 3-d can 
one speak sensibly about a well-defined net accreted angular momentum per unit rest-mass, 
jin- On the basis of these contrasts, we conclude that axisymmetric simulations, while useful 
for preliminary investigations of various initial configurations, are of only limited utility in 
investigating the long term behavior of MHD turbulent disks. 



3.2. The coronal mechanism 



The dominant mechanism for bringing flux to the event horizon in simulation VDOmSd is 
one that operates primarily in the corona, and not in the disk accretion flow proper. Although 
most of the flux is brought in during the first 5000M in time, there is a slow increase in 
flux throughout the remainder of the simulation. The majority of this flux appears to be 
transported through the continuing operation of the coronal mechanism. The key aspect is 
that magnetic field is carried inward by infalling low density fluid. Although this coronal 
mechanism operates outside the disk where most of the mass accretion happens, it, too, 
depends on angular momentum transport to proceed, and the basic physical mechanism for 
angular momentum transport is likewise the same: magnetic torques. Wherever the vertical 
field is perturbed so as to create a radial component, orbital shear acting on the radial field 
in turn creates toroidal field. Moreover, whenever the orbital frequency decreases outward, 
the relative signs of the radial and toroidal fields are always such as to make the r-0 element 
of the magnetic stress tensor, b^b^/ATT, negative. Such a stress transports angular momentum 
outward, tending to accentuate inflow, and therefore growth in the radial field component; 
this is, of course, the mechanism of the magneto-rotational instability. 
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An example of the coronal mechanism in action is given in Figure [8] In this figure, 
color contours show the strength of the azimuthally averaged Maxwell stress per unit mass, 
overlaid with field lines derived from A^{r,6). For clarity we focus on those field lines 
that still pass through the equator in the accretion disk. The first frame, corresponding to 
t = 6000M, features a field that runs from the equator vertically upward through the disk 
and then down toward the hole before reversing and moving radially outward. The strong 
stress per unit mass in the vicinity of the hairpin carries angular momentum away from the 
associated matter, enabling it to move inward. In the next frame at t = 6200, the hairpin 
moves in radially toward the black hole. The next two frames show a close up view of the 
innermost part of the accretion fiow at times t = 6440 and 6460M. By the time of the 
first one, the inner portion of the hairpin has reached and crossed the horizon. Only 20M 
later, the inside hairpin has reconnected, forming a closed field loop. Thus, the net result 
of the coronal mechanism at this stage is to attach those fieldlines forming the exterior of 
the hairpin to the black hole. These fieldlines then extend out into the corona, ultimately 
attaching to the disk at much larger radius. 

Although net local infiow of magnetic field is a necessary condition for fiux accretion, 
it is not a sufficient one; the net fiux changes only when reconnection occurs. Local motion 
alone is insufficient because magnetic fiux is an essentially global concept, so local motion is 
not enough to determine net fiux transport. This distinction is emphasized by contrasting 
the local magnetic infiow rate with global measures. We define the local rate by 

V^{r,9) = ^ [ d<i)V'B' (17) 
ivr J 

averaged over the portion of the simulation during which the trapped fiux grows slowly 
(5000-2 X lO^M). This is shown in Figure [oj for which the radius was chosen to be 15M. As 
this figure shows, the local fiux infiow speed is considerably greater just outside the funnel 
(Itt — 6'| ~ 0.27r), where it is ~ 0.5-1 x 10~^, than in the midplane, where it is nearly two 
orders of magnitude smaller. 

This local rate of fiux infiow is, however, much larger than the global rate of fiux accu- 
mulation. The global rate is dA^/dt, which can be inferred from Figure |4j For the purpose 
of this figure, the magnetic fiux was separated into the fiux piercing the midplane, $, and the 
fiux through the upper hemisphere of the black hole, \E'. During the early stages of the simu- 
lation (before 5000M), d-^/dt ~ 1 x 10~^ but this rate of increase diminishes at later times 
by at least an order of magnitude. Thus, the local infiow rate, V$(6' ~ 0.27r, O.Svr) ~ 1 x 10~^ 
at r = 15M during the latter 3/4 of the simulation, is as much as two orders of magnitude 
greater than the net global rate. In other words, local field infiow alone does not suffice to 
create global fiux infiow. 



- 16 - 



Likewise, reconnection on its own is also a necessary condition for flux accretion, but 
not a sufficient one. In order for flux accretion to occur, reconnection must result in a 
global rearrangement of the field topology. This fact is illustrated by the field structures 
displayed in Figure [8j In that figure, a narrow current sheet lies between the two sides of 
the hairpin, permitting reconnection to happen there. In order to transfer magnetic flux 
inward through the accretion flow and increase the net flux in the funnel, reconnection must 
preferentially destroy the side of the magnetic hairpin that lies closest to the midplane. This 
is possible because reconnection can also occur across the equator, something that the funnel 



field cannot do. Figure [TO] illustrates the idea. Coronal hairpin field loops bring field down 
to the inner disk, returning along the disk back out to large radius where they eventually 
pass through the equator and connect to another hairpin in the opposite hemisphere. The 
innermost bend of these hairpins can attach to the horizon. When this happens, the field 
appears to be almost a closed loop, with one end of the loop inside the horizon (e.g., the 
two inner-disk field lines in the gray region of the diagram). If the two oppositely-directed 
portions of one of these loops can meet at intermediate radius, the equatorial flux jumps 
inward. This inward jump drastically shortens the time until the entire field loop can accrete 
onto the black hole. When that happens, flux of one sign is removed from the horizon above 
the equator and flux of the other sign is removed from below the equator, leaving a net 
increase in flux piercing each horizon hemisphere. 

An example of a reconnection event within VDOmSd that leads to this sort of global 



flux rearrangement is shown in Figure 11 In this figure, color contours denote the radial 
component of the CT-field, indicating the direction of the fieldline, which are again 
overlaid with white contours derived from the azimuthally averaged A^{r,9). Initially (t = 
14600M), fieldlines connect large radii to the black hole event horizon through the funnel 
region (deep blue region below midplane) and then proceed back out into the disk close to 
the midplane (yellow region below midplane). These fieldlines cross the midplane at large 
radius and proceed back to the black hole north of the midplane (light blue region north 
of midplane), before finally proceeding out to large radius through the funnel region (red 
region north of midplane). At time t = 14620M, the oppositely directly fieldlines that lie 
just above and below the midplane have reconnected at r ~ 6M, and the resulting closed 
fieldloop (which is still attached to the black hole) is rapidly accreted. This process transfers 
flux from the outer disk to r ~ 6M (where the reconnection event occurs) and then down to 
the black hole as the field loop is accreted, arriving at the black hole with the outer edge of 
the loop. 



The coronal mechanism resembles the flux advection mechanism proposed by Rothstein 



& Lovelace (2008) in the sense that most flux transport takes place outside of the disk body. 



and that flux in the corona moves inward because magnetic stresses send angular momentum 
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off to infinity. They emphasize the contrast between the turbulent disk and non-turbulent 
flow outside of the disk. We likewise observe that large-scale, relatively laminar flows within 
the low-density corona easily move fleld about; flux-freezing is an excellent approximation 
there. On the other hand, there are also several points of contrast between our results and 
the model of Rothstein & Lovelace (2008). They argue that in the disk corona, deflned as 
the region where /3 ~ 1, the component of the magnetic stress tensor could transport 
angular momentum outward as effectively as the component to which attention is more 
often directed, allowing the disk surface to move inward. Although this is similar to what 
we see in our simulation during the initial avalanche phase, it is different from the dynamics 
prevailing at later times. During the accretion phase of the simulation, the coronal fleld 
lines develop MRI-like bends at high altitudes above the disk (>> H), so that the inflow 



occurs in the corona rather than along the disk surface. In addition, as shown by Figure 12 



which contrasts the azimuthal averages of —b^b^/p and —b^b^/p, there is very little net 
angular momentum flux per unit mass in the polar angle direction except in the funnel (at 
latitudes not too far from the midplane, b^b^p ~ b^b^). Although the magnitude of is often 
comparable to T^, it has very nearly the same probability of having either sign, so that its 
mean value is small. On the funnel edge, is more likely to take angular momentum away 
from the disk but has magnitude on average somewhat smaller than T^, while in the outer 
corona, on average it brings angular momentum toward the midplane. 



3.3. Inflow in the disk body 

The majority of the net flux on the black hole is brought there by the coronal mechanism 
during the initial 5000M of time. However, the flux on the horizon does increase over the last 
1.5 X lO^M and, while the coronal mechanism continues to operate during this interval, it is 
possible that part of the flux might be brought in through the accretion disk itself. As the 
right panel of Figure |9] makes plain, the local magnetic inflow rate V$ is roughly two orders of 
magnitude smaller in the midplane than it is at the funnel edge. Nonetheless, its magnitude 
is still consistent with the total flux inflow rate during the latter part of the simulation. It 
may therefore be possible for flux advection directly associated with the accretion flow in 
the disk body to account for some part of the flux accumulation during the quasi-steady 
accretion phase. In this section we attempt to estimate the size of that contribution. 

We begin by comparing the rate of magnetic flux accretion to the rate of mass accretion. 
If magnetic flux and mass moved inward in lock-step, the functions 



M{r,t) = [ dt'M{r,t')/M^ 
Jo 



(18) 
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and 

Air, t) = A^{r, t, 7r/2)M<^( w, 0, tt/2) (19) 
would be identical. These quantities are normalized to the initial mass, Mq, and the initial 



flux; Tmax IS the outer radius of the initial torus. In Figure [T3| we show the time-dependence 
of both these quantities at two fiducial radii, just outside the event horizon (left panel) and 
at r = 20M (right panel). The former, of course, tells about the flux attached directly to the 
black hole and the mass deposited there. The difference between the two plots is the mass 
and flux within the accretion disk between the horizon and r = 20M; here we will refer to 
this region as the "inner disk." 



In evaluating Figure 13, one must allow for an important distinction between the his- 
tories of mass and magnetic flux: the mass of the black hole must increase monotonically 
(and the mass of the inner disk, although not required to do so, generally does), but the 
net magnetic flux can increase or decrease, both by magnetic reconnection and by the radial 
motion of closed loops. Thus, the inflow of magnetic flux is highly intermittent, and the time- 
derivative of flux within some radius can easily be negative. On short timescales, the ratio 
of the rate of flux accretion to mass accretion can vary by at least an order of magnitude, 
as well as fluctuate in sign. For this reason, a sensible comparison of the rate of magnetic 
flux inflow relative to that of mass can be made only in regard to long-term time averages. 
In addition, in order to evaluate the efficiency of this process in real disks, we must restrict 
these averages to times when the inner disk was in approximate inflow equilibrium: for this 
simulation, that means the latter half of the data, from t = lO'^M to t = 2 x lO'^M. Within 
this span, we smooth A{t) by taking a running average lOOOM wide and then compose the 
ratio of the change in ^ to 7W from lO^M to 2 x lO^M using the smoothed data. We find 
that AA/AAi ~ 0.47 at the horizon, and ~ 0.50 at r = 20M. In other words, relative to 
the initial amount embedded in the accretion flow, magnetic flux moves inward roughly half 
as fast as mass, and these rates are the same at both the horizon and r = 20M. If the mass 
is in a state of inflow equilibrium, the equality of these two ratios shows that the magnetic 
flux in the inner disk is also in equilibrium. 

This ratio of inflow rates may be contrasted with the local ratio of magnetic flux to 



mass. Looking beyond the short-timescale fluctuations in the plots of A{t) in Figure 13 
one sees that the average values of the flux at the horizon and the flux contained within 
r = 20M are very similar. In other words, the vertical magnetic flux through the inner 
disk, $(r = 20M), contributes very little to A. Time-averaging over the last lO^M of data, 
we find that it is only ~ 0.86% of the initial flux, and it appears that what flux there is is 
confined to the outer part of the inner disk. By contrast, averaged over the same period, 5% 
of the initial mass can be found in the inner disk. Thus, in the inner disk, the time-averaged 
flux/mass ratio is only 17% of the initial value. 



- 19 - 



To reconcile this flux/mass ratio with the accreted flux/mass ratio, there are only two 
possibilities: either the magnetic flux moves in at the same rate (or slower) than the mass 
and (at least) 2/3 of the flux reaching the horizon between lO^M and 2 x lO^M arrived 
there via a different route (most likely the coronal mechanism), or the magnetic flux on 
average moved inward three times faster than the mass. We view the former alternative as 
much more plausible. Note, however, any flux delivery from outside the inner disk must be 
balanced by flux losses from the inner disk to the horizon because the magnetic inflow rate 
at r = 20M is nearly the same as at the horizon. 

Next we further explore the nature of the magnetic flux equilibrium in the inner disk. 



using two different views of this quantity. Figure 14 shows the distribution of the vertical flux 
function $(r) time-averaged over the last 5000M of the simulation and normalized to the 
initial flux value. For comparison, we also plot in that figure both the initial flux distribution 
and the final mass distribution in the disk, both likewise normalized to the initial total. The 
radial derivative d^/dr shows the magnitude and sign of the vertical magnetic field piercing 
the equator. We see that, like the mass, vertical flux has spread away from its initial location, 
both in and out, but, as already mentioned, a much smaller part of the net magnetic flux 
than the mass resides in the inner disk. 

The time-dependence of the magnetic flux in the inner disk over the last lO^M in time is 



illustrated in Figure 15 As this figure shows, not only is the magnetic flux contained within 
the inner disk a rather small part of the total, it fluctuates in time, frequently going negative. 
Little sign of any long-term trend can be seen, consistent with our contention that the inner 
disk magnetic flux has reached a statistical equilibrium. The frequent sign changes of the 
local net flux suggest that the poloidal magnetic field lines in the inner disk predominantly 
close within the inner disk. 

Moreover, the magnetic field corresponding to the net flux is a very small fraction of the 



total field within the disk. Figure [16j shows the azimuthally averaged value of B through 
the equator at the end time as a function of radius, along with the initial value. The MRI 
has generated considerable field of both signs throughout the disk. Local regions of the disk 
have a net vertical flux, but their contrasting signs lead to little contribution to the total 
disk flux. The global net flux is present only as a slight positive excess moving both inward 
and outward with time. 

To reach the state described by these figures, most of the large-scale magnetic connec- 
tions to the disk matter must have been reconnected away. But this should come as no 
surprise because we have already seen that a significant part of the original flux has moved 
hy reconnection from the disk to the horizon. The matter in the inner disk has little net 
flux because it has been largely bypassed through the coronal mechanism. Thus, while the 
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MRI has created a large effective turbulent viscosity in the sense that considerable mass has 
accreted from the location of the initial torus, very little net flux has moved with it. In 
other words, within the disk, turbulence created by the MRI transports angular momentum 
rapidly and drives accretion, but very little flux moves with the accreting mass. 



3.4. Strength of the Funnel Field 



Lastly, we raise the question of whether in the very long-term the magnetic flux on the 
horizon would continue to grow. If it is regulated by a combination of the magnetic and gas 
pressure near r = 6M , and these are in an approximately time-steady state, one might expect 



that further flux accumulation on the horizon would have to stop. Figure [17] supports this 
idea. This figure shows that the magnetic pressure inside the funnel, the magnetic pressure 
in the inner accretion flow, and the gas pressure in the inner accretion flow are well correlated 
with one another in a temporal sense, even while all three change by more than two orders 
of magnitude over the course of the simulation. The time-dependent data show fluctuations 
that are so large, however, that the much slower accumulation of flux seen in Figure [13] is 



well below the noise, so our data on magnetic and gas pressures cannot answer any questions 
about long-term saturation of magnetic flux attached to the horizon. 



That the magnetic and gas pressures should be closely related is an old idea (Rees 



et al. 1982 MacDonald & Thorne 1982), but there is also a long history of controversy 



about whether the magnetic field pressure in the funnel should be closer in magnitude to the 



(poloidal) magnetic pressure or the gas pressure in the inner disk. For example, Ghosh & 



Abramowicz (1997) and Livio et al. (1999) argued that as the poloidal magnetic field strength 



in the disk is subthermal, the funnel field strength should be regulated by it rather than the 
gas pressure. In the simulation presented here, the total, i.e. poloidal plus toroidal magnetic 
pressure in the disk is smaller than the gas pressure for the duration of the simulation (i.e. 
the field remains subthermal), and the magnetic pressure from the field in the funnel lies 
between these two candidate regulators, sometimes closer to one and sometimes closer to the 
other. Thus, the funnel field is consistently stronger than the poloidal field in the disk. 

The relevance of the inner disk poloidal field to the strength of the funnel field seems 
to be rather limited based upon what we observe about the process carrying flux to the 
horizon. The fleld in the funnel is nearly radial, and its total intensity is determined by the 
poloidal flux that is delivered to the black hole over the course of the simulation. This flux 
delivery system is the coronal mechanism, and so it should not be surprising that the poloidal 
flux within the accretion disk itself plays at most a secondary role (e.g., by determining the 
rate of reconnection and accretion of flux loops). Instead, simple pressure- matching in the 
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inner regions of the accretion flow — be it gas, (total, mostly toroidal) magnetic or radiation 
pressure — appears to be important in regulating how much flux may be attached to the black 
hole. We speculate that when the flux on the horizon has a total field intensity larger than 
would be consistent with the gas and magnetic pressure in the nearby accretion flow, the 
rate of reconnection along the funnel fleld boundary rises, so that the field approaching the 
black hole is entirely in closed loops and does not add to the net fiux on the horizon. 

4. Summciry, Discussion &c Conclusions 

In this paper we investigate the evolution of an accretion torus embedded in a large-scale 
vertical magnetic field orbiting a Schwarzschild black hole, with a view toward studying how 
magnetic fiux moves relative to the accretion fiow. The simulation stretched over 2.0 x lO^M 
in time, corresponding to 80^2^^ at the initial torus pressure maximum, long enough to 
establish inflow equilibrium in the inner disk for the second half of the simulation. Of 
particular interest is how the net vertical fleld evolves, and whether or not a fleld distribution 
consistent with the formation of jets or winds can develop. We trace the evolution of the net 
poloidal flux distribution with a particular focus on net flux that becomes attached to the 
black hole. Our primary result is that a signiflcant fraction of the initial flux — ~ 27% — is 
brought to the black hole horizon, even though only a rather smaller fraction of the initial 
mass — ~ 19% is accreted. The flux attached to the horizon supports a coherent poloidal 
field within the evacuated axial funnel, a requirement for the creation of a Blandford-Znajek 
type jet (if the black hole rotates, which it did not in this simulation). The mass and flux 
are carried inward through distinct mechanisms: the mass by turbulent stresses within the 
accretion disk, and the flux by large-scale motions in the low density corona. 

Most of the global magnetic flux motion is mediated by a novel mechanism that we 
have dubbed the coronal mechanism. Rather than the gradual "diffusive" process that had 
been the focus of most previous discussions of magnetic flux motion in accretion flows, 
in this mechanism the flux is brought directly down to the horizon as the net flux jumps 
discontinuously when large-scale magnetic loops reconnect across the equator. The same 
orbital shear that creates angular momentum transport in the disk body by correlating radial 
and azimuthal magnetic field also acts in the corona; the difference is that in the corona the 
resulting stress leads to the formation of large-scale loops stretched rapidly inward rather 
than to MRl-driven turbulence. These stretched loops can reconnect to oppositely-directed 
field at much smaller radius than their footpoint, leading to sudden macroscopic radial 
changes in the location of magnetic flux. Although the coronal mechanism acts particularly 
rapidly during the initial transient phase of the simulation, it continues to dominate flux 
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motion throughout the simulation, including the long period during which the inner disk is 
in approximate inflow equilibrium. We might therefore reasonably expect it to be a property 
of actual accretion disks. 

The same reconnection that creates global relocation of net flux simultaneously creates 
closed fleld loops within the disk body. For this reason, the flux/mass ratio within the 
accretion flow is suppressed by a factor of order unity relative to what it is in the initial 
state. Earlier global disk simulations have, in most cases, assumed zero net flux for reasons 
of simplicity and computational convenience; the coronal mechanism might, in fact, make 
this assumption a reasonable approximation to the magnetic fleld in typical disks. The 
suppression of flux/mass in the disk body also makes conventional advection of flux in 
direction association with mass accretion relatively inefficient. 



As seen in previous simulations (particularly Beckwith et al. 2008a), if the overall fleld 
topology at the black hole is dipolar, the funnel fleld can be relatively immune to reconnection 
and hence long-lived. We conjecture that the funnel fleld amplitude is determined by pressure 
balance with the gas and (total, mostly toroidal) magnetic pressures in the inner disk and 
coronal regions. Future experiments could test this hypothesis by running much longer 
simulations with more available net fleld to see if the fleld strength levels off or continues to 
build. 

As with any simulation, our results are subject to some uncertainty due to the lim- 
itations inherent in a numerical solution. One concern is the importance of reconnection 
to the overall flux evolution. In our numerical simulation, reconnection occurs at the grid 
scale when oppositely signed fleld components are brought together. Although leaving the 
microscopic rate of such an important process to numerical effects is a concern, the overall 
motions that lead to the formation of current sheets and subsequent reconnection area driven 
by events on much larger scales. In this sense, we might argue that the rate of reconnection 
seen in the simulation is relatively independent of our gridscale. We have not, however, 
demonstrated that our results are numerically converged, although the global fleld motions 
in the axisymmetric simulations carried out at two different resolutions were very similar. 
Furthermore, the thermodynamics of these simulations are incomplete because energy lost 
to numerical magnetic reconnection and numerical cancellation of fluid velocities is not cap- 
tured; conversely, we do not account for radiation in any way, either as a cooling agent or 
as a contributor to the pressure. One thing that is clear, regardless of any other numerical 
limitation, is that axisymmetric simulations are of limited utility. In addition to the con- 
straints associated with the anti-dynamo theorem, the vertical fleld channel modes remain 
dominant throughout the simulation, giving a qualitatively distinct evolution at late time. 
Three dimensional simulations are essential for studying long term, steady-state behavior. 
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To place our somewhat surprising results in context, it is worthwhile contrasting them 
with previous suggestions about how magnetic flux moves through accretion flows. These 
are quite disparate, for they have in general been based on approximate analytic arguments. 
We examine three principal approaches, of which one has received much more attention than 
the others. 

The greatest amount of attention has been given to a picture in which the matter moves 
inward by an unspecified angular momentum transport process (called "viscosity", but not 
thought to be literally that) while the magnetic field diffuses relative to the matter through 
another mechanism called "turbulent resistivity", but similarly unattached to specific mi- 
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Livio 2001). The ratio between this effective viscosity and resistivity (the nominal Prandtl 
number), would then determine whether the flux mostly moves inward with the mass, or 
instead diffuses outward relative to the mass fast enough to avoid much net inward motion. 

This approach rests upon two little-examined assumptions: First, that the overall flow 
can be regarded as if it were laminar and time-steady, and second, that the behavior of 
the underlying turbulence can be reduced to parameters (an advection rate and a diffusion, 
i.e., resistivity, coefficient) whose values are independent of the magnetic field. With these 
assumptions, field motion in the disk body governs field motion far away from the disk, 
and one may describe the field evolution in the language of linear diffusion. Unfortunately, 
neither assumption applies to real accretion flows. 

It is now well-established that angular momentum transport in accretion disks is due to 
turbulent magnetic stresses driven by the MRI. Consequently, the magnetic field structure 
is neither laminar nor time-steady. Even in the corona, where the MRI per se does not 
occur, the same basic physics leads (as we have earlier discussed) to irregular field motions. 
There is therefore no direct connection between local field motion deep inside the disk (e.g., 
resistive diffusion) and the position of the corresponding field lines far away. Moreover, as 
we have emphasized, flux is a global concept, not a local one; the coronal mechanism, for 
example, could never be described by a local theory of this sort. 

The second assumption, that the evolution of the field can be described in terms of 
simple resistive and viscous diffusion operating on a background field gradient is similarly 
problematic. Accretion is driven by nonlinear magnetic stresses; because the effective inflow 
velocity arises only from an average of a strongly fluctuating velocity field, it is non-local in 
both time and space. Similarly, the breakdown in ideal MHD happens primarily by driven 
reconnection, the rate of which depends on the structure and magnitude of both the magnetic 
and velocity fields. As a result, neither the mean rate at which magnetic field is carried from 
place to place nor any tendency for its structure to smooth is either independent of or linear 
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in the field strength. 

In fact, the formulation in terms of a competition between "viscous" and "resistive" 
diffusion breaks down in an even more fundamental way. Where there is net flux passing 
through the disk body, the very same MHD turbulence that transports angular momentum 
simultaneously decouples the flux from the matter through turbulence-driven reconnection. 
Field lines associated with flux passing through the equator develop bends that readily break 
off as closed loops. The matter moves inward with the closed loops while the flux stays in 
place. Thus, in dramatic contrast to the conventional prediction, rapid "viscous" accretion 
can, and apparently does, co-exist with largely stationary flux, while the "resistivity" that 
permits this decoupling has no particular impact on smoothing the field distribution. 

Moreover, there is relatively little net flux passing through the disk because reconnection 
associated with the coronal mechanism detaches it efficiently. In other words, most of the 
flux motion occurs in an "end-run" that bypasses the disk altogether, so that it has rather 
little to do with any activity in the disk body, "viscous" , "resistive" , or other. 

Local smoothing of field structures can occur, but it is questionable whether one can de- 



fine an overall diffusion coefficient to describe its rate. For example, Guan & Gammie (2009) 
have recently attempted to quantify magnetic diffusion in the context of MRI-driven MHD 
turbulence using shearing box simulations. They imposed a sinusoidally varying vertical 
field with an amplitude above the background due to the MRI-turbulent flow and observed 
the subsequent decrease in this mode's Fourier power. The peak vertical field amplitudes 
studied vary from 10-80% of the initial toroidal field strength. They found that the decay 
time for the imposed feature is several tens of Q~^, with the decay rate an increasing func- 
tion of the perturbation amplitude. These results indicate that local field gradients can be 
smoothed within the turbulent flow, but it is unclear to what degree the inferred average 
diffusion coefficient depends on the specific situation studied. In particular, in accretion 
disks, the net field, the part determining the flux, is, especially in the disk body, small in 
magnitude compared to the fluctuating turbulent field. When the dynamics are strongly 
nonlinear, spreading of structures in this small net field may proceed in a way very different 
from spreading of larger amplitude structures; in fact, the motion of the net field may have 
much more to do with the dynamics of the fluctuating field than to its own disposition. 

For all these reasons, therefore, we see no useful way to describe the results of the present 
simulation in the traditional language of arbitrarily specified transport coefficients applied 
to gradients of the time-averaged magnetic field. 



Two other concepts have recently been studied. One of these (Rothstein & Lovelace 



2008 Lovelace et al. 2009) in some respects is a blend of the "advection/diffusion" picture 
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and the coronal mechanism. On the one hand, it uses the conventional formalism of fixed 
transport coefficients; on the other hand, it relies heavily on field line motions in the corona. 
These two approaches are united by supposing that the MRI and its resulting turbulence 
are suppressed in the corona, so that it suffers neither "turbulent viscosity" nor "turbulent 
resistivity". Consequently, coronal motion is dominated by vertical transport of angular 
momentum, in a magneto-centrifugal wind when the effective Prandtl number is large, and 
in a Poynting flux-dominated outflow otherwise. The result is that (for most parameters), 
the upper surface of the disk moves inward, carrying the footpoints of the large-scale flux 
with it, while the flow in the midplane is outward. 

As we have previously discussed, the fundamental physical element of the MRI, the 
creation of substantial magnetic stresses by orbital shear, does operate within the corona, and 
strongly — this is the foundation of the coronal mechanism. Despite these coronal motions, 
this is not a significant source of angular momentum transport for the main accretion disk. 
Although the instantaneous amplitude of the magnetic stress is typically sizable, outside 
the funnel its sign fluctuates so that we see no significant net vertical transport of angular 
momentum from the disk. In addition, the net fluid radial velocity in the midplane, although 
small compared to its rms value, is inward. Nonetheless, our results are consistent with their 
suggestion that magnetic stresses carrying angular momentum outward through the corona 
can be effective in moving field lines inward. Because the mass density in the corona is so low, 
removing a comparatively small amount of angular momentum can lead to large-scale field 
motions. Because the subsequent motions are not dominated by turbulence, flux freezing 
remains effective within the coronal plasma. 



In the second of the two recent suggestions for how flux moves in accretion flows. Spruit 



& Uzdensky (2005) propose that the flux motion is not controlled by a simple gradient, 
but rather by the dynamics of intermittent field bundles. These authors suggest that mag- 
netized patches could accrete relatively rapidly from angular momentum losses through a 
wind. Further, if such fields were strong enough to suppress the MRI, then turbulence, 
and subsequently magnetic reconnection, would be greatly reduced for those patches. While 
such magnetized patches would not themselves lead to a significant net magnetic flux pass- 
ing through the disk per se, they could accumulate at the central black hole until the field 
strength itself suppressed further accretion. Such a scenario is similar to that of the mag- 



netically arrested accretion model of Narayan et al. (2003). Strong nonaxisymmetric field 



structures appear in a simulation of Igumenshchev (2008) where a strong field is contin- 



ually injected at the outer boundary, providing some support for the [Spruit Uzdensky 



(2005) concept. In our simulation, we do see nonaxisymmetric variations in the vertical 
flux throughout the disk, but not in the form of local pockets of intense vertical flux, nor 
do we see either any significant disk wind or any interruption of the accretion flow due to 



"magnetic arrest." Nonetheless, we agree with Spruit & Uzdensky (2005) in emphasizing 
the intermittency of the flow and the potential importance of nonaxisymmetries. 



The study carried out here is, of course, only a first step, and the results presented here 
point to several avenues of investigation with future numerical experiments. For example, 
simulations that continue for longer times and begin with net field at larger radii can better 
explore the efficacy of the coronal mechanism. In particular longer simulations with more 
available net field are required to see if the field strength in the funnel levels off or continues 
to build over time. Simulations in the Kerr metric with non-zero spin parameters could 
also study directly how magnetic flux evolution and the coronal mechanism relate to jet 
launching and coUimation. The requirements of three spatial dimensions and increased 
spatial resolution will continue to make these investigations challenging. 
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Fig. 1. — Initial torus and field configuration. White contours denote magnetic field lines, 
color contours density distribution (left panel) and the gas (3 parameter (right panel). 
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Fig. 2. — Initial evolution of the flow in simulation VDOmSd showing the advection of the 
vertical field from the initial location of the torus to the black hole. White contours are 
magnetic field lines, color contours show the density distribution. From left to right, the 
panels correspond to times t = 1000, 1500, 2000, 2500M. The color scale is the same as in 
Figure [1} By t = 2500M, the funnel region has become strongly magnetic pressure dominated 
and two large scale channel modes have formed within the initial torus. 
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Fig. 3. — Growth and subsequent break up of the vertical field MRI channel mode in sim- 
ulation VDOmSd. White contours are magnetic field lines, color contours show the density 
distribution. From left to right, the panels correspond to times t = 3000, 3500, 4000, 4500M. 
The color scale is the same as in Figure [TJ By t = 4500M, the disk has become turbulent; 
large scale MRI modes persist within the corona. 
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Fig. 4. — Total flux integrated from the axis along the horizon to the equator and then 
out along the equator to the outer boundary. The equatorial component, and the horizon 
component, are shown separately. The vertical scale is in units of the total initial flux. 
The total flux is unchanged through the evolution; losses in the equatorial flux are accounted 
for in net flux through the horizon. 
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Fig. 5. — Azimuthally- averaged gas density (left panel) and plasma f3 (right panel) at time 
2 X lO^M in simulation VD0m3d. 
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Fig. 6. — Contours of A^{r, 9), displayed both in color and as poloidal field lines (isosurfaces 
of A^{r, 6)) at time 2 x lO^M in simulation VDOmSd. 
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Fig. 7. — Comparison between results from the two-dimensional (VD0m2d) and three- 
dimensional (VD0m3d) simulation. Top panel shows the flux through the horizon, ^ for 
VD0m2d (solid hne) and VDOmSd (dashed hne) . The middle panel shows the net mass ac- 
cretion rate through the horizon in units of the initial torus mass. The bottom panel shows 
the specific angular momentum carried into the hole as a function of time for simulations 
VD0m2d and VDOmSd. The strong coherent radial fields in the two-dimensional system 
create large fluctations in these quantities. 
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Fig. 8. — Evolution of a magnetic "hairpin" in the corona during the turbulent steady state. 
Color contours show the azimuthal average of stress per unit mass, —b^b^/p, and are overlaid 
with the poloidal field structure (white lines), which are kept fixed between the panels. From 
left to right, the top panels show times t = 6000, and 6200M (top row). The bottom row 
shows the near-hole region at times t = 6440 and 6460M. During this time a reconnection 
event takes place within the field hairpin. 



-37- 




2^10 I I I I I I I I I I I I I I I 

0.46 0.48 0.50 0.52 0.54 



Fig. 9.— Azimuthal mean of -V''B^{15M), averaged over the time span 5000-1.5 x lO^M at 
r = 15M. Two 9 ranges and scales are used: a scale appropriate for the corona and excluding 
only the funnel region (top), and only the disk proper plotted on a scale appropriate to the 
value there. 
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Fig. 10. — Schematic of the field fine structure in the coronal mechanism. Field lines within 
the coronal are carried in toward the black hole, forming a hairpin-like structure. When the 
hairpin connects to the horizon fiux is added to the funnel field and opposite-signed fiux is 
added to the disk (shaded region). Reconnection across the equator (dashed line) allows this 
field to form loops that accrete. Accretion of those field loops results in an increase of the 
net horizon fiux. 
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Fig. 11. — Reconnection event that occurs in the disk and leads to the rapid transport of 
poloidal magnetic flux from the outer disk to the black hole. Color contours show and 
are overlaid with the poloidal field structure (white lines), which are kept fixed between the 
panels. From left to right the panels show t = 14600, 14620, and 14640M. 
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Fig. 12. — Azimuthal average of —h^h^/p (left panel) and —b^b^/p (right panel) at the same 
time as shown in Fig. [sj Regions with —b^b^/p > 1 are colored dark red, regions with 
—b^b^/p < 10"^ (including —b^bif,/p < 0) are deep blue. The sign of —b^b^/p is adjusted so 
that positive values correspond to electromagnetic angular momentum flux directed away 
from the disk. The azimuthal average of —b^b^ is much smaller than that of —b^bfp everywhere 
outside the funnel. 
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Fig. 13. — Fractional accretion of magnetic flux A (solid curve) and mass Ai (dashed curve) 
as functions of time. Because we stored shell-integrated mass accretion rates every IM in 
time, but the full 3-d data required to compute the vector potential only every lOM, the 
latter curve has greater time resolution than the former. (Top) At r — 2.1M. (Bottom) At 
r = 20M. 
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Fig. 14. — The net vertical flux through the equator $(r), integrated from the black hole 
horizon outward, averaged over the last 5000M of the simulation and normalized to the inital 
total value ("Vertical Flux, $"). The initial $ distribution is also shown (dashed line). The 
net flux has spread through the disk over the course of the simulation. For comparison, the 
mass distribution M(r) is shown (labeled "Mass") for the end time, t = 2 x 10^ M, hkewise 
normalized to its initial total value. 
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Fig. 15. — Space-time contours of the net vertical flux through the equator $(r), integrated 
from the black hole horizon outward for the inner disk during the second half of the sim- 
ulation, normalized to the total initial flux. Never very large in magnitude, the integrated 
magnetic flux in this region fluctuates in sign. 
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Fig. 16. — Azimuthal average of the absolute value of \B^\ (solid line) and the average of 
(dot-dashed hne) along the equatorial plane at i = 2 x 10^ M (solid hne) along with 
the initial distribution (dashed line). The vertical field remains relatively strong but 
experiences fluctuations in sign throughout the disk body. The net flux is a small positive 
bias in the vertical field value. 
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Fig. 17. — Azimuthally- averaged magnetic pressure as a function of time at two locations: 
inside the funnel near its base (r = 6M, 9 = 0.0447r; solid curve); and in the inner accretion 
flow (r = 6M, 6 = 7r/2; dotted curve), and azimuthally-averaged gas pressure in the inner 
accretion flow (also r = 6M, = ■n/2] dashed curve). 



